home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / MATH1 / ERF4.PAS < prev    next >
Pascal/Delphi Source File  |  1985-04-03  |  2KB  |  81 lines
  1. program erfd4;        { -> 334 }
  2.  
  3. { evaluation of the gaussian error function }
  4.  
  5. var    x,er,ec        : real;
  6.     done        : boolean;
  7.  
  8. external procedure cls;
  9.  
  10. function erf(x: real): real;
  11. { infinite series expansion of the Gaussian error function }
  12.  
  13. const    sqrtpi        = 1.7724538;
  14.     t2        = 0.66666667;
  15.     t3        = 0.66666667;
  16.     t4        = 0.07619048;
  17.     t5        = 0.01693122;
  18.     t6        = 3.078403E-3;
  19.     t7        = 4.736005E-4;
  20.     t8        = 6.314673E-5;
  21.     t9        = 7.429027E-6;
  22.     t10        = 7.820028E-7;
  23.     t11        = 7.447646E-8;
  24.     t12        = 6.476214E-9;
  25.  
  26. var    x2,sum        : real;
  27.     i        : integer;
  28.  
  29. begin
  30.   x2:=x*x;
  31.   sum:=t5+x2*(t6+x2*(t7+x2*(t8+x2*(t9+x2*(t10+x2*(t11+x2*t12))))));
  32.   erf:=2.0*exp(-x2)/sqrtpi*(x*(1+x2*(t2+x2*(t3+x2*(t4+x2*sum)))))
  33. end;    { function erf }
  34.  
  35. function erfc(x: real): real;
  36. { complement of error function }
  37. const    sqrtpi        = 1.7724538;
  38.  
  39. var    x2,v,sum    : real;
  40.  
  41. begin
  42.   x2:=x*x;
  43.   v:=1.0/(2.0*x2);
  44.   sum:=v/(1+8*v/(1+9*v/(1+10*v/(1+11*v/(1+12*v)))));
  45.   sum:=v/(1+3*v/(1+4*v/(1+5*v/(1+6*v/(1+7*sum)))));
  46.   erfc:=1.0/(exp(x2)*x*sqrtpi*(1+v/(1+2*sum)))
  47. end;        { function ercf }
  48.  
  49. begin        { main }
  50.   cls;
  51.   done:=false;
  52.   writeln;
  53.   repeat
  54.     write('Arg? ');
  55.     readln(x);
  56.     if x<0.0 then done:=true
  57.     else
  58.       begin
  59.     if x=0.0 then
  60.       begin
  61.         er:=0.0;
  62.         ec:=1.0
  63.       end
  64.     else
  65.       begin
  66.         if x<1.5 then
  67.           begin
  68.         er:=erf(x);
  69.         ec:=1.0-er
  70.           end
  71.         else
  72.           begin
  73.         ec:=erfc(x);
  74.         er:=1.0-ec
  75.         end    { if }
  76.     end;
  77.     writeln('X= ',x:6:2,', Erf= ',er:7:4,', Erfc= ',ec)
  78.       end    { if }
  79.     until done
  80. end.
  81.